# ------------------------------------- #
# Replication code for:
#
# Rathbun, Brian, Christopher Sebastian Parker, and Caleb Pomeroy "Separate but Unequal: Ethnocentrism and Racialization 
# Explain the 'Democratic' Peace in Public Opinion," American Political Science Review.
#
# This script reproduces the Qualtrics survey results presented in the main text and appendices.
# ------------------------------------- #

# --- libraries --- #
library(psych)
library(ggeffects)
library(ggplot2)
library(texreg)
library(ltm)
library(plyr)
library(dplyr)
library(reshape2)
library(estimatr)
set.seed(1912)

# --- set working directory
setwd("~/Downloads/replication_files/")

# --- load data
qualtrics_survey <- readRDS("data/qualtrics_survey.rds")
nrow(qualtrics_survey) # N = 1626, as mentioned in the paper's text 

# --- ethnocentrism measure (from Neuliep and McCroskey)
ethno_items <- cbind(qualtrics_survey$ethno1,
                     qualtrics_survey$ethno2,
                     qualtrics_survey$ethno3)
cronbach.alpha(ethno_items) # alpha = 0.68, as footnoted in main text
ethno_fit <- fa(ethno_items, nfactors=1, rotate="varimax", scores=TRUE, fm="pa")
ethno_fit$loadings # SS loadings 1.32, as footnoted in main text
#qualtrics_survey$ethno_nm_factor <- ethno_fit$scores[,1]

# --- prep data
qualtrics_survey$Y <- qualtrics_survey$nukes_strike
qualtrics_survey$W <- ifelse(qualtrics_survey$nukes_regime == "dem", 1, 0)
qualtrics_survey$Z_0 <- ifelse(qualtrics_survey$nukes_race == "only", 1, 0)
qualtrics_survey$Z_1 <- ifelse(qualtrics_survey$nukes_race == "white", 1, 0)
qualtrics_survey$Z_2 <- ifelse(qualtrics_survey$nukes_race == "nonwhite", 1, 0)
qualtrics_survey$E <- ifelse(qualtrics_survey$ethno_nm_factor > median(qualtrics_survey$ethno_nm_factor), 1, 0)
qualtrics_survey$MI <- ifelse(qualtrics_survey$mil_int > median(na.omit(qualtrics_survey$mil_int)), 1, 0)
qualtrics_survey$SDO <- ifelse(qualtrics_survey$sdo > median(na.omit(qualtrics_survey$sdo)), 1, 0)
qualtrics_survey$RWA <- ifelse(qualtrics_survey$rwa > median(na.omit(qualtrics_survey$rwa)), 1, 0)
qualtrics_survey$CON <- ifelse(qualtrics_survey$ideology > median(na.omit(qualtrics_survey$ideology)), 1, 0)
qualtrics_survey$R <- ifelse(qualtrics_survey$race_resent_factor > median(na.omit(qualtrics_survey$race_resent_factor)), 1, 0)

# --- figure A6
{plot_demos <- qualtrics_survey
  plot_demos$ids <- rownames(plot_demos)
  plot_demos$age_plot <- ifelse(qualtrics_survey$age < 30, "18-29",
                                ifelse(qualtrics_survey$age >= 30 & qualtrics_survey$age < 40, "30-39",
                                       ifelse(qualtrics_survey$age >= 40 & qualtrics_survey$age < 50, "40-49",
                                              ifelse(qualtrics_survey$age >= 50 & qualtrics_survey$age < 60, "50-59",
                                                     ifelse(qualtrics_survey$age >= 60 & qualtrics_survey$age < 70, "60-69", "70+")))))
  plot_demos <- reshape2::melt(plot_demos[,c("ids", "gender", "age_plot", "ideology")], id.vars = "ids")
  plot_demos$variable <- plyr::revalue(plot_demos$variable,
                                       c("gender" = "Gender",
                                         "age_plot" = "Age",
                                         "ideology" = "Political Ideology"))
  qualtrics_survey$ids <- rownames(qualtrics_survey)
  ed_df <- data.frame(ids = qualtrics_survey$ids,
                      variable = "Education",
                      value = qualtrics_survey$education)
  p1 <- 
    ggplot(subset(plot_demos, variable == "Gender"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Female", "Male")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p2 <- 
    ggplot(subset(plot_demos, variable == "Age"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p3 <- 
    ggplot(ed_df, aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("College\nGraduate", "Less than\nCollege")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p4 <- 
    ggplot(subset(plot_demos, variable == "Political Ideology"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Extremely\nLiberal", "2", "3", "4", "5", "6", "Extremely\nConservative")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  
  gridExtra::grid.arrange(p1,p2,p3,p4, ncol = 1)
}


# --- figure A7
{p1 <- 
  ggplot(qualtrics_survey, aes(ethno_nm_factor)) +
  geom_histogram(bins = 20, color = "black", fill ="grey90") +
  labs(x = "Ethnocentrism Factor Scores\n(Neuliep and McCroskey)", y = "Count", tag = "A") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 15, margin = margin(t=15)),
        axis.title.y = element_text(size = 14, margin = margin(r=15)),
        panel.spacing.x = unit(1.5, "lines"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"),
        plot.tag = element_text(face = "bold", size = 22),
        plot.tag.position = c(.06, .95))
p2 <- 
  ggplot(qualtrics_survey, aes(ethno_nm_additive)) +
  geom_bar(color = "black", fill ="grey90") +
  labs(x = "Ethnocentrism Additive Scores\n(Neuliep and McCroskey)", y = "Count", tag = "B") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 15, margin = margin(t=15)),
        axis.title.y = element_text(size = 14, margin = margin(r=15)),
        panel.spacing.x = unit(1.5, "lines"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"),
        plot.tag = element_text(face = "bold", size = 22),
        plot.tag.position = c(.06, .95))
gridExtra::grid.arrange(p1, p2, nrow = 1)}


# --- qualtrics analysis: full sample ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W, data = qualtrics_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95 <- lm_robust(Y ~ W * Z_0, data = qualtrics_survey, 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0, data = qualtrics_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A4
screenreg(list(lm_dem_95, lm_white_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W, data = qualtrics_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90 <- lm_robust(Y ~ W * Z_0, data = qualtrics_survey, 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0, data = qualtrics_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the full sample at alpha .05 and .10 levels; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_white_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_white_95$coefficients[4], lm_nonwhite_95$coefficients[4]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_white_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_white_95$conf.high[4], lm_nonwhite_95$conf.high[4]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_white_95$conf.low[2], lm_nonwhite_95$conf.low[2],  lm_white_95$conf.low[4], lm_nonwhite_95$conf.low[4]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_white_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_white_90$conf.high[4], lm_nonwhite_90$conf.high[4]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_white_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_white_90$conf.low[4], lm_nonwhite_90$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "full_sample")

# --- qualtrics analysis: high ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

# --- table A5
screenreg(list(lm_dem_95_he, lm_white_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .05)
summary(lm_white_90_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from high ethno respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_white_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_white_95_he$coefficients[4], lm_nonwhite_95_he$coefficients[4]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_white_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2], lm_white_95_he$conf.high[4], lm_nonwhite_95_he$conf.high[4]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_white_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_white_95_he$conf.low[4], lm_nonwhite_95_he$conf.low[4]),
             ci_90_high = c(lm_dem_90_he$conf.high[2], lm_white_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[2], lm_white_90_he$conf.high[4], lm_nonwhite_90_he$conf.high[4]),
             ci_90_low = c(lm_dem_90_he$conf.low[2], lm_white_90_he$conf.low[2], lm_nonwhite_90_he$conf.low[2], lm_white_90_he$conf.low[4], lm_nonwhite_90_he$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- vary the ethnocentrism threshold
# in the main text, we present results with ethnocentrism binned at the median
# here, we simply note that this result is robust beyond this median split
# nonwhite information significantly eliminates the democratic peace effect for approximately 57% of the sample, as mentioned in the paper's text

screenreg(lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, ethno_nm_factor >= quantile(qualtrics_survey$ethno_nm_factor, c(.43))), 
                    subset = (Z_0 == 1 | Z_2 == 1),
                    alpha = .05))

# --- qualtrics analysis: low ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A6
screenreg(list(lm_dem_95_le, lm_white_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from low ethno respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_white_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_white_95_le$coefficients[4], lm_nonwhite_95_le$coefficients[4]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_white_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2], lm_white_95_le$conf.high[4], lm_nonwhite_95_le$conf.high[4]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_white_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_white_95_le$conf.low[4], lm_nonwhite_95_le$conf.low[4]),
             ci_90_high = c(lm_dem_90_le$conf.high[2], lm_white_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[2], lm_white_90_le$conf.high[4], lm_nonwhite_90_le$conf.high[4]),
             ci_90_low = c(lm_dem_90_le$conf.low[2], lm_white_90_le$conf.low[2], lm_nonwhite_90_le$conf.low[2], lm_white_90_le$conf.low[4], lm_nonwhite_90_le$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Qualtrics Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure 2
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type, scales = "free_x") +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}

# --- qualtrics analysis: compare ethnocentrism to other individual differences ---#
# here, we consider various individual differences to assess whether the ethnocentrism results are unique 
# --- median split for differences in ATE
summary(mod_ethno <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
summary(mod_mi <- lm_robust(Y ~ W, data = subset(qualtrics_survey, MI == 1), 
                            subset = Z_0 == 1,
                            alpha = .05))
summary(mod_sdo <- lm_robust(Y ~ W, data = subset(qualtrics_survey, SDO == 1), 
                             subset = Z_0 == 1,
                             alpha = .05))
summary(mod_rwa <- lm_robust(Y ~ W, data = subset(qualtrics_survey, RWA == 1), 
                             subset = Z_0 == 1,
                             alpha = .05))
summary(mod_con <- lm_robust(Y ~ W, data = subset(qualtrics_survey, CON == 1), 
                             subset = Z_0 == 1,
                             alpha = .05))
summary(mod_rr <- lm_robust(Y ~ W, data = subset(qualtrics_survey, R == 1), 
                            subset = Z_0 == 1,
                            alpha = .05))
# --- table A10
screenreg(list(mod_ethno, mod_mi, mod_sdo, mod_rwa, mod_con, mod_rr))

# --- interaction between individuals differences and the eliminated effect
# here, the bottom-most coefficient in each table is the eliminated effect x individual diff quantity of interest
# the p-value associated with that coefficient is in the "Pr(>|t|)" column, as reported in the appendix text.
summary(mod_ethno <- lm_robust(Y ~ W * Z_0 * E, data = qualtrics_survey, 
                               subset = (Z_0 == 1 | Z_2 == 1),
                               alpha = .05))
summary(mod_mi <- lm_robust(Y ~ W * Z_0 * MI, data = qualtrics_survey, 
                            subset = (Z_0 == 1 | Z_2 == 1),
                            alpha = .05))
summary(mod_sdo <- lm_robust(Y ~ W * Z_0 * SDO, data = qualtrics_survey, 
                             subset = (Z_0 == 1 | Z_2 == 1),
                             alpha = .05))
summary(mod_rwa <- lm_robust(Y ~ W * Z_0 * RWA, data = qualtrics_survey, 
                             subset = (Z_0 == 1 | Z_2 == 1),
                             alpha = .05))
summary(mod_con <- lm_robust(Y ~ W * Z_0 * CON, data = qualtrics_survey, 
                             subset = (Z_0 == 1 | Z_2 == 1),
                             alpha = .05))
summary(mod_rr <- lm_robust(Y ~ W * Z_0 * R, data = qualtrics_survey, 
                            subset = (Z_0 == 1 | Z_2 == 1),
                            alpha = .05))

# --- qualtrics analysis: full sample (with covariates) ---#
# here, we simply re-run the above models with the addition of right hand side covariates
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95 <- lm_robust(Y ~ W * Z_0  + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A7
screenreg(list(lm_dem_95, lm_white_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W  + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90 <- lm_robust(Y ~ W * Z_0  + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0  + gender + race + pid + ideology + education + age, data = qualtrics_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the full sample at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_white_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_white_95$coefficients[length(lm_white_95$coefficients)], lm_nonwhite_95$coefficients[length(lm_nonwhite_95$coefficients)]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_white_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_white_95$conf.high[length(lm_white_95$coefficients)], lm_nonwhite_95$conf.high[length(lm_nonwhite_95$coefficients)]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_white_95$conf.low[2], lm_nonwhite_95$conf.low[2],  lm_white_95$conf.low[length(lm_white_95$coefficients)], lm_nonwhite_95$conf.low[length(lm_nonwhite_95$coefficients)]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_white_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_white_90$conf.high[length(lm_white_90$coefficients)], lm_nonwhite_90$conf.high[length(lm_nonwhite_90$coefficients)]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_white_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_white_90$conf.low[length(lm_white_90$coefficients)], lm_nonwhite_90$conf.low[length(lm_nonwhite_90$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "full_sample")

# --- qualtrics analysis: high ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W  + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_he <- lm_robust(Y ~ W * Z_0  + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0  + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A8
screenreg(list(lm_dem_95_he, lm_white_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90_he <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the high ethno respondents at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_white_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_white_95_he$coefficients[length(lm_white_95_he$coefficients)], lm_nonwhite_95_he$coefficients[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_white_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2], lm_white_95_he$conf.high[length(lm_white_95_he$coefficients)], lm_nonwhite_95_he$conf.high[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_white_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_white_95_he$conf.low[length(lm_white_95_he$coefficients)], lm_nonwhite_95_he$conf.low[length(lm_nonwhite_95_he$coefficients)]),
             ci_90_high = c(lm_dem_90_he$conf.high[2], lm_white_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[2], lm_white_90_he$conf.high[length(lm_white_90_he$coefficients)], lm_nonwhite_90_he$conf.high[length(lm_nonwhite_90_he$coefficients)]),
             ci_90_low = c(lm_dem_90_he$conf.low[2], lm_white_90_he$conf.low[2], lm_nonwhite_90_he$conf.low[2], lm_white_90_he$conf.low[length(lm_white_90_he$coefficients)], lm_nonwhite_90_he$conf.low[length(lm_nonwhite_90_he$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- qualtrics analysis: low ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A9
screenreg(list(lm_dem_95_le, lm_white_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(qualtrics_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the low ethno respondents at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_white_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_white_95_le$coefficients[length(lm_white_95_le$coefficients)], lm_nonwhite_95_le$coefficients[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_white_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2], lm_white_95_le$conf.high[length(lm_white_95_le$coefficients)], lm_nonwhite_95_le$conf.high[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_white_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_white_95_le$conf.low[length(lm_white_95_le$coefficients)], lm_nonwhite_95_le$conf.low[length(lm_nonwhite_95_le$coefficients)]),
             ci_90_high = c(lm_dem_90_le$conf.high[2], lm_white_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[2], lm_white_90_le$conf.high[length(lm_white_90_le$coefficients)], lm_nonwhite_90_le$conf.high[length(lm_nonwhite_90_le$coefficients)]),
             ci_90_low = c(lm_dem_90_le$conf.low[2], lm_white_90_le$conf.low[2], lm_nonwhite_90_le$conf.low[2], lm_white_90_le$conf.low[length(lm_white_90_le$coefficients)], lm_nonwhite_90_le$conf.low[length(lm_nonwhite_90_le$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results w/ covariates modeled
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Qualtrics Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure A8
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type, scales = "free_x") +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}


# --- qualtrics analysis: main effect of race ---#
# as noted in the dataverse appendix, we find little evidence of a main effect of race on support for a strike
t.test(subset(qualtrics_survey, nukes_race == "nonwhite")$Y,
       subset(qualtrics_survey, nukes_race == "white")$Y)


# --- qualtrics analysis: nonwhite identifying respondents ---#
# here, we re-estimate the primary experimental results for the nonwhite identifying respondents in our sample
# however, note: the experiment was not designed with sufficient statistical power to identify these effects.
# --- all nonwhite identifying respondents
nrow(subset(qualtrics_survey, race == "nonwhite"))
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W, data = subset(qualtrics_survey, race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95 <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B1
screenreg(list(lm_dem_95, lm_white_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W, data = subset(qualtrics_survey, race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90 <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_white_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_white_95$coefficients[4], lm_nonwhite_95$coefficients[4]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_white_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_white_95$conf.high[4], lm_nonwhite_95$conf.high[4]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_white_95$conf.low[2], lm_nonwhite_95$conf.low[2],  lm_white_95$conf.low[4], lm_nonwhite_95$conf.low[4]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_white_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_white_90$conf.high[4], lm_nonwhite_90$conf.high[4]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_white_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_white_90$conf.low[4], lm_nonwhite_90$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "full_sample")


# --- nonwhite identifying respondents high in ethnocentrism
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B2
screenreg(list(lm_dem_95_he, lm_white_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 1 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents above ethno median, at alpha .05 and .10 levels; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_white_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_white_95_he$coefficients[4], lm_nonwhite_95_he$coefficients[4]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_white_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2], lm_white_95_he$conf.high[4], lm_nonwhite_95_he$conf.high[4]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_white_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_white_95_he$conf.low[4], lm_nonwhite_95_he$conf.low[4]),
             ci_90_high = c(lm_dem_90_he$conf.high[2], lm_white_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[2], lm_white_90_he$conf.high[4], lm_nonwhite_90_he$conf.high[4]),
             ci_90_low = c(lm_dem_90_he$conf.low[2], lm_white_90_he$conf.low[2], lm_nonwhite_90_he$conf.low[2], lm_white_90_he$conf.low[4], lm_nonwhite_90_he$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- nonwhite identifying respondents low in ethnocentrism
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- white mediator arm (alpha = .05)
summary(lm_white_95_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B3
screenreg(list(lm_dem_95_le, lm_white_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- white mediator arm (alpha = .10)
summary(lm_white_90_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                                 subset = (Z_0 == 1 | Z_1 == 1),
                                 alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0, data = subset(qualtrics_survey, E == 0 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents below ethno median, at alpha .05 and .10 levels; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE", "ACDE", "Eliminated Effect", "Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_white_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_white_95_le$coefficients[4], lm_nonwhite_95_le$coefficients[4]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_white_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2], lm_white_95_le$conf.high[4], lm_nonwhite_95_le$conf.high[4]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_white_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_white_95_le$conf.low[4], lm_nonwhite_95_le$conf.low[4]),
             ci_90_high = c(lm_dem_90_le$conf.high[2], lm_white_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[2], lm_white_90_le$conf.high[4], lm_nonwhite_90_le$conf.high[4]),
             ci_90_low = c(lm_dem_90_le$conf.low[2], lm_white_90_le$conf.low[2], lm_nonwhite_90_le$conf.low[2], lm_white_90_le$conf.low[4], lm_nonwhite_90_le$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = White", "Country Race = Nonwhite", "Country Race = White", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Qualtrics Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure B1
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type, scales = "free_x") +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}

